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Exploring the hydrodynamic boundary of a surface by approaching a colloidal sphere and measuring 
the occurring drag force is a common experimental technique. However, numerous parameters like the 
wettability and surface roughness influence the result. In experiments these cannot be separated easily. 
For a deeper understanding of such surface effects a tool is required that predicts the influence of different 
surface properties. In this paper we present computer simulations based on the lattice Boltzmann method 
of a sphere submerged in a Newtonian liquid. We show that our method is able to reproduce the theoretical 
predictions for flat and noninteracting surfaces. In order to provide high precision simulation results the 
influence of finite size effects has to be well controlled. Therefore we study the influence of the required 
system size and resolution of the sphere and demonstrate that even moderate computing resources allow 
the error to be kept below 1%. 
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1. Introduction 

The miniaturization of technical devices down to submicrometric sizes has led to the development of so- 
called microelectro-mechanical systems (MEM S) whi ch are now commonly applied for chemical, bio- 
logical and technical applications (see Tabeling (2005)). A wide variety of microfiuidic systems was de- 
veloped including gas chromatography systems, electrophoretic separation systems, micromixers, DNA 
amplifiers, and chemical reactors. Further, microfiuidic experiments were used to answer fundamental 
questions in physics including the beh avior of single m o lecules or pa r ticles in fluid flow or the validity 
of the no-slip boundary condition (see iTabelind (12005b : lLauga et al.1 (12005b ). A violation of the latter 
in sub-micron sized geometri es was found in very well con t rolled experiments in recent yea r s. Sinc e 




then, mostly experimental dLauga e t al. (2005); Crai g et al.l (12001 ): iTrethewav and Meinhart (2004a); 
Cheng and Giordand(l2002T);IChoi et al.ld2003h:lBaudrv and Charlaixld200lh : ICottin-Bizonne et al.ld2002l) 



Vinogradova and Yakubov (2003); Churaev et al. ( 1984)), but also theoretical studies like lVinogr adova 



Gennes ( 20021); Bazant and Vinogradova ( 2008 ) , as well as computer simula tions (see Sued 
Barr at and BocqueJUii^ Ciep lak et al.l d200lb : iThompson and Troianl dl997h ) have been per- 



formed to improve our understanding of slippage. The topic is of fundamental interest because it has 
practical consequences in the physical and engineering sciences as well as for medical and industrial 
applications. Interestingly, also for gas flows, often a slip length much larger than expected from classi- 
cal theory can be observed. Ex t ensive revie ws of the slip phenomenon have recently been published by 
Lauga et all d2005b . lNeto et al.l d2005b and iBocquet and Barratl d2007l) . 

Boundary slip is typically quantified by the slip length b. This concept was proposed by iNaviei 



( 1 18231) . He introduced a boundary condition where the fluid velocity at a surface is proportional to the 
shear rate at the surface (at x — xq), i.e. 



v z {x ) 



dv z (x) 



(1.1) 



In other words, the slip length b can be defined as the distance from the surface where the relative flow 
velocity vanishes. 

The experimental investigation of slip can be based on different setups. Early reports of slippage 
measured the mass flow Q through a pipe with radius R with a controlled pressure gradient and compare 
it to the theoretical values for Poiseuille flow. The slip length b can then be derived from the ratio 
between the measured and the theoretical mass flow as 



Qp = u 4b . 

^Poiseuille ^ 



(1.2) 



Schnelll dl956h used such a method to measure the slip length in coated glass capillaries and f ound 
b = 1 — lO^m. O ther groups found only a slip length of b = 10 — 70nm (see lLauga et al.1 (12005b and 
Neto et al. (2005)). The problem is that in such experiments a deviation of the capillary radius R leads 
automatically to a measured slip. With other words one can not distinguish between a shift in the 
boundary position or an actual slip phenomenon. 

Other experimental methods follow the flow field by introducing tracer particles into the flow. The 
basic assumption off all these methods is that the tracers do have the same velocity as the flow, and 
t hat they do no t distu rb the flow. An example for such a method is the double focus cross correlation 
( iLumma et alj d2003l) ). By having two laser beams at a fixed distance it is possible to record every 
particle that crosses one of the focal points of the beams. The technique rests in the premise that only 
a small number of labeled particles are simultaneously located in an effective focal volume of the order 
of 10~ 15 1. Therefore the time cross-correlation can be used to determine the average time a particle 
needs to cross the second focus after it crosses the first one. Since the focus of the laser beams can be 



locate d very precisely one has an accurate measurement of the flow velocity at a given point. iPit et al 



(2000) 1 lave applied the cross correlation method for hexadecane flowing along a hydrocarbon lyophobic 
smooth surface and found a slip length of 400nm. IVinogradova et al.1 (|2009) refined the method and 
determined a slip length for water and NaCl aqueous solutions of less than lOOnm for a hydrophobic 
polymer channel which is independent on the shear rate and the salt-concentration. For hydrophilic 
surfaces no measurable slip was detected. 

A second example of such a particle based method is micro particle image velocimetry (micro PIV). 
The method is very similar to the cross correlation method, but instead of just taking into account single 
events of a particle crossing the focus of a laser beam, one takes a series of pictures and correlates 
them to each other. For macroscopic phenomena this technique is typically realized by a CCD video 
camera but on the microscopic level, more sophisticated methods are needed. This includes the use 
of two lasers with different color to illuminate the tracers. Fu rther, the optical setup is c rucial , so tha t 



the focal plane of the optical set up is in the right position. Tret hewav and Mei nhart (2002, 2004b) 
applied m icro PIV for water in hydr ophobic glass channels and measured slip lengths of up to 1/im. 
However I Joseph and Tabelingl ( 20051) found slip lengths for water of less than lOOnm, stating that this 
is the minimal resolution of the method, i.e., that it is doubtful whether there is any slip. 

Both methods give the average velocity distribution of the tracer particles in the flow field. This 
measured flow profile can then be fitted by a theoretical flow profile. For a Poiseuille flow, driven by the 



pressure gradient with the viscosity \i the channel width d and the slip length b on both sides of the 
channel, the flow profile in x direction reads as 
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1 dp 
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(1.3) 



As with the simple mass flow measurements reported on above a principle problem of these tracer based 
methods is that it is not possible to distinguish between a shift in the boundary position and intrinsic 
slip, i.e. in Eq. 1 1 .31 the slip length b and the (effective) channel width d cannot be decoupled without 
a variation of d. Another drawback is in the basic assumption of the tracer models, namely that the 
flow velocity is equal and undisturbed by the particles. However on the micro scale the particle-particle 
interaction and particle-wall interaction might become an issue. Van der Waals and electrostatic forces 
between the fluid or surface and the particle can lead to a significant change of the particle trajectory. 

Another class of experimental methods to determine slip is based on the measurement of the drag 
between two surfaces or the force on a particle moving towards a boundary. Very popular is the 
modification of an atomic force microscope (AFM) by adding a silicon sphere to the tip of the cantilever. 
A sketch of such an experiment is shown in Fig. [T] While moving the surface towards the sphere, the 
drag force F can be measured with a high precision. Since the typical distance d between the sphere 
and the surface as well as the approaching velocity v are small the simple Reynolds theory for the 
lubrication force 

F Re = 6nn\R 2 /d, (1.4) 

can be applied, where jj. is the dynamic viscosity of the fluid and R the radius of the sphere. Eq. II. 41 is 
valid for a no-slip sphere approaching a no-slip surface. To measure the amount of slip at the wall, b, the 
drag force is compared with its the oretical value for a slip surface, as reported bv lVinogradova and Yakubov 
d2003l) : I Vinogradoval (119961 1 1995b . A correction /* is applied to Eq. 11 .41 that takes into account the sur- 
face properties 

F b = fF{d). (1.5) 

In case of a sur face with the slip le ngth b and a vanishing slip on the surface of the sphere, the correction 
/* is given by IVinogradoval (11995b as 
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There are some possibly problematic limitations of this setup. Those include the hydrodynamic influ- 
ence of the cantilever, a twist of the cantilever or uncontrolled roughness on the sphere or the approached 
surface. Further, the exact velocity of the sphere is hard to control since the drag force and the bending 
of the cantilever can lead to an acceleration of the sphere. This contributes to a deviation from the ideal 
case and causes sophisticated corrections of the approaching velocity to be required. However, if the 
experimental setup is well controlled, it allows a very precise measurement and allows one to distinguish 
between a shift in the boundary position or slip. 



In general it should be noted that a significant dispersion of the slip me asurements rep < 
Tature can be observed even in similar experimental systems (see ILauga et al.l q2005); 



orted in the 



Neto etal 



literature 

d2005h ). This is due to the large amount of unknown or uncontroll ed parameters . For example, ob- 
served slip lengths vary between a few nanometres as reported by IChuraev et al. I (11984 and microm- 
eters as shown in iTrefhewav and M einhart (2004a) and while some authors like IChoi et al. (2003); 




FIG. 1, A sketch of a slip measurement based on a modified atomic force microscope (AFM). For technical reasons the surface 
is moved with the velocity v towards the sphere which is attached to the cantilever of the AFM. The laser is used to measure the 
bending of the cantilever which can be related to the drag force F. 



Zhu and Granick d200ll): Craig et al. (2001 find a dependence of the slip on the flow velocity, others 
like lCheng and Giordanol(l2002h : lTretheway and Meinhartl d2004al) do not. 

The large variety of different experimental results to some extent has its origin in surface-fluid 
interactions. Their properties and thus their influence on the experimental results are often unknown 
and difficult to quantify. In addition there are many influences that lead to the same effect in a given 
experimental setup as intrinsic boundary slip - as for example a fluid layer with lower viscosity than the 
bulk viscosity near the boundary. Unless one is able to resolve the properties of this boundary layer it 
cannot be distinguished from true or intrinsic slip. Such effects can be categorized as apparent slip. 

In the literature a large variety of different effects leading to an apparent slip can be found. How- 
ever, detailed experimental studies are often difficult or even impossible. Here computer simulations 
can be utilized to predict the influence of effects like surface wettability or roughness. Most recent 
simulations apply molecular d ynamics and report increasing slip with decreasing liquid density (see 
Thompson and Robbinsl ( 1990l) ) or liquid-solid interactions (see Cieplak et alj( 2001 ); Nagavama and C heng 
(2004)), while slip decreases with increasing pressure (see Barrat and Bocquet ( 1999)). These simula- 
tions are usually limited to some tens of thousands of particles, length scales of the order of nanometres 
and time scales of the orde r of nanoseconds. A lso, shear rates are usually orders of magnitude higher 
than in any experiment (see lLauga et alJ (120051) ). Due to the small accessible time and length scales of 
molecular dynamics simulations, mesoscopic simulation methods like the lattice Boltzmann method are 
highly applicable for the simulation of microfluidic experiments. For a simple flow setup l ike Po i seuille 
or Couette flow several i nvestigations of slip models have been published by iNie et alJ ( 120021) : ISucci 
d2002h : ISbragaglia etaD J2006I) : lHarting et ail (12006b : iHyvaluoma and Hartingl d2008l) . but investiga- 
tions for more complex flows are rare. 

In this paper lattice Boltzmann simulations of AFM based experiments are presented. We focus on 
demonstrating that our method is able to reproduce the theoretical prediction for a simple no-slip case 
and investigate its limits. It is then possible to use the results to investigate the influence of different 
parameters. Such parameters can include an interaction between the sphe re and the boundary , surface 
roughness or hydrophobicity. While the case of roughness is investigated in lKunert et al. (2010), further 
studies will be shown in future publications. Microscopic origins of slip caused by the molecular details 
of the fluid and surface cannot be treated with a continuum method and are beyond the scope of the 
current paper. 



The remainder of this paper is arranged as follows: after this introduction we describe the theoret- 
ical background and our simulation method. Then we show that the method is able to reproduce the 
theoretical predictions with great accuracy. Further, we investigate the influence of finite size effects. 
We determine the minimum system size and resolution of the discretization required to push the finite 
size effects below an acceptable limit and demonstrate the limits of the simulation method. 



2. Theoretical background 

In this section the commonly used theory is summarized. The Reynolds lubrication force is given by 
Eq |1.4| For larger distances d between the surface of the sphere and the approached boundary this 
force does not converge towards the Stokes drag force Fg t = 6npL\R for a sphere moving freely in a 
fluid. Therefore this simple equation fails in cases of larger separations d, where the Stokes f orce is 



not su fficie ntly sma l l to be neglected. The system can be described accurately by the theory of Maude 



dl96lh and lBrennen (1196 II) . The base of the theory is a solution for two spheres approaching each other 
with the same rate. By transforming the coordinates, applying symmetry arguments and setting the 
radius of one of the spheres to infinity one arrives at a fast converging sum for the drag force acting on 
the sphere: 

FMa = 67E/iV/tti, (2.1) 



with 
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where cosh^ = (d — R) /R. The term given by X\ cannot be treated analytically. Thus, we evaluate X\ 
numeric ally with a con vergence of 10~ 10 . A more practical approximation of ( 12.11 ) is given in the same 
paper bv lMaudddl96lh : 

, x f9R \ 

F{d)=6%ilRv\-- + l\ (2.2) 

Here one can easily see that the force converges towards the Stokes force for an infinite distance d and 
towards the Reynolds lubrication ( 1 1 .4b for small separations d. 

As stated above, to measure the slip length b experiments apply a correction /* that takes into 
account the surface properties so that fj, = f*F{d). The correction function /* is given by Eq. 11.61 
This equation is valid for a perfectly flat surface with finite slip, but it does not allow one to distinguish 
between slip and other effects like surface roughness. Therefore it is of importance to perform computer 
simulations which have the advantage that all relevant parameters can be changed independently without 
modifying anything else in the setup. Thus, the influence of every single modification can be studied in 
order to present estimates of the influence on the measured slip lengths. The first step in this process 
is to validate the simulation method and to understand its merits and flaws. In general most computer 
simulations suffer from the fact that only a small system can be described and that one is usually not able 
to simulate the whole experimental system in full detail. For this reason it is mandatory to understand 
which resolution of the problem is required to keep finite size effects under control and to cover the 
important physics correctly. 




In Fig.|2]a sketch of the simulated system is shown. A sphere embedded in a surrounding fluid is simu- 
lated. While the sphere approaches a surface the force acting on it is recorded. We perform simulations 
with different sphere radius R and different system length L to investigate finite size effects. For sim- 
plicity we set the x and z dimension to the same value L and keep the propagation dimension y constant 
at 512 lattice units. A typical approach to limit finite size effects is to use peri odic boundary con ditions. 
However in such a system the sphere then interacts with its periodic images. lHasimotoldl959h gives a 



theoretical solution for the drag force of a sphere in a periodic array as it appears if all boundaries are 
periodic: 
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(2.3) 



Here, a = R/L is the ratio between the radius of the sphere R and the system length L. The case 
we are dealing with in this paper is different and more complex due to a broken symmetry caused by 
the approached surface. For the approximation given by (12.3b the main contribution of the periodic 
interaction is between the periodic images in front and behind the sphere. Due to the rigid boundary, 
these are not present in our case. 

Besides a finite simulated volume most simulation methods utilize a finite discretization of the sim- 
ulated objects, i.e. the sphere in our case. This means that the finite size and the resolution influence the 
result of a simulation. However, it is usually possible to limit the influence of finite size effects and the 
loss of accuracy by discretization if those errors are known and taken into account properly. Therefore, 
we study the finite size effects in a simulation of a sphere in a periodic system approaching a rigid no 
slip boundary and investigate how different resolutions of the sphere influence the force acting on it. 



3. Simulation Method 

The simulation method used to study microfluidic devices has to be chosen carefully. While Navier- 
Stokes solvers are able to cover most problems in fluid dynamics, they lack the possibility to include the 
influence of molecular interactions as needed to model boundary slip. Molecular dynamics simulations 
(MD) are the best choice to simulate the fluid-wall interaction, but computer power today is not sufficient 
to simulate length and time scales necessary to achieve orders of magnitude which are relevant for 
experiments. However, boundary slip with a slip length b of the order of man y molecular diameters 



a has been studied with molecular dynamics simulations by various authors like lThompson and Troian 



dl997l) : ICieplak et all d200ll) : ICottin-Bizonne et all d2004 : lBaudrv and Charlaixl d200lb . 

In this paper we use the lattice Boltzmann method, where one discretizes the Boltzmann kinetic 
equation 

d F 
uV*+-' 



dt 



«(x,u,f ) = Q 



(3.1) 



on a lattice. n(x,u,f) indicates the probability to find a single particle with mass m and velocity u at the 
time t and position x. F accounts for external forces. The derivatives represent simple propagation of a 
single particle in real and velocity space whereas the collision operator £2 takes into account molecular 
collisions in which a particle changes its momentum due to a collision with another particle. Further, 
the collision operator drives the distribution n towards an equilibrium distribution n eq . 

In the lattice Boltzmann method time, positions, and velocity space are discretized on a lattice in the 
following way. The distribution n is only present on lattice nodes x^. The velocity space is discretized so 
that in one discrete timestep 8t the particles travel with the discrete velocities c, towards the nearest and 
next nearest neighbours x# + c,5f . Since a large proportion of the distribution stays at the same lattice 
node a rest velocity Co is required. Co represents particles not moving to a neighboring site. In short 
we operate on a three dimensional grid with 19 velocities (i = 0..18) which is commonly referred to as 
D3Q19. After the streaming of the population density, the population on each lattice node is relaxed 
towards an equilibrium such that mass and momentum are conserved. It can be shown b y a Chapman - 
Enskog procedure that such a simulation method reproduces the Navier-Stokes equation. ISuccil (12001b . 
In the lattice Boltzmann method the time t is discretized in time steps dt, the position x is discretized in 
units of distance between neighbouring lattice cells, and the velocity u is discretized using the velocity 
vectors c,-. These form the natural lattice units of the method which are used in this paper if not stated 
otherwise. 



The implementation we are using originates from I Verberg and L add (2001). It applies a so called 
multi relaxation time collision operator. The distribution n is transformed via a transformation matrix T 
into the space of the moments rrij of the distribution. A very accessible feature of this approach is that 
some of the moments mj have a physical meaning as for example the density 



mo 



(3.2) 



or the momentum in each direction x,y,z, 



2,3 =ei,2,3E n ' C '' 



(3.3) 



ei 2 3 is the unit vector in Cartesian directions. The moments m; relax with an individual rate Sj towards 
the equilibrium mj 1 . The equilibrium distribution m eq conserves mass mo and momentum nt\ ,2,3 and 
is a discretized version of the Maxwell distribution. Thus the lattice Boltzmann equation ( 13. U can be 
written as 



n(xjfc + c,-5f,Cj,f + St) - n(x k ,d,t) 
= T i Sj[mj(x k ,t)-mf). 



(3.4) 



The multi relaxation time approach has several advantages compared to other lattice Boltzmann schemes. 
These include a higher precision at solid boundaries and the direct accessibility of the moments which 
represent physical properties. The latter can be util ized to easily implement thermal fluctuations or 
external and internal forces (see lDiinweg et al.1 d2007D ). 



A feature of the implementation we are using is the possibility to si mulate particles suspended in 



fluid. T he simulation method is d e scribed extensively in t he literature (seelR. Verberg and A. J. C. Ladd 



(200d); IVerberg and Ladd! d200ll) : lDunweg et all (2007): Komnik et al] (12004). Therefore only a brief 



description is given here. The movement of the particles is described by a simple molecular dynamics 
algorithm. However it should be noted here that we simulate a sphere moving with a constant velocity. 
Therefore any forces acting on the sphere do not influence its movement. The fluid-particle interaction 
is achieved by the solid-fluid boundary interaction acting on the surface of the sphere. When the sphere 
is discretized on the lattice all lattice sites inside the sphere are marked as boundary nodes with a moving 
wall boundary condition. This boundary condition at the solid-fluid interface is constructed in such a 
way that there is as much momentum transferred to the fluid as required for the fluid velocity to match 
the boundary velocity v& of the particle. The center of mass velocity of the particle and the rotation are 
taken into account. This way the transfered momentum and thus the hydrodynamic force acting on the 
sphere are known. Technically speaking a link bounce back boundary condition is implemented for the 
solid nodes together with a momentum transfer term. The link bounce back implies that the distributions 
that would move inside the boundary with the velocity c, are reversed in direction with opposite velocity 

Cjfc. 

n(x + Cj8t,d,t + 8t) = nhL + Ci8t,Ck,t) + ^ a ^ VfeC ' (3.5) 

, c» x /x 2a Ci p\bCi 
n(x,c k ,t + 8t) = «(x,c,-,f)H ~ (3.6) 

Here, a Ci are weight factors taking into account the different lengths of the lattice vectors. 

While the center of mass of the sphere moves, new lattice nodes become part of the particle, while 
others become fluid. Therefore particles do not perform a continuous movement but rather small jumps. 
After each jump the fluid is out of equilibrium but relaxes very quickly back to the quasi static state. In 
order to average out statistical fluctuations imposed by the discrete movement of the particle, one has to 
average the recorded force over several time steps. We choose to average over intervals of 999 steps. 

If not stated otherwise the simulation parameters are v = 0.001, jx = 0.1 and the radius is varied 
between R = 4 and R = 16. The approached boundary is a plain no-slip wall which is realized by a mid 
grid bounce back boundary condition. 

Along the open sides periodic boundary conditions are applied so that the sphere can interact with 
its mirror leading to the to be avoided finite size effects. As noted in Fig.[2]the length of the system in 
x and y direction is L. The size of the simulation volume is varied to explore the influence of finite size 
effects. 



4. Results 

In this contribution we vary the system length L and the radius of the sphere R. A major contribution to 
the finite size effects is the interaction of the sphere with its periodic image. Therefore a larger system 
length should reduce this effect dramatically. However when the hydrodynamic influence of the wall 
becomes larger finite size effects become smaller. This can be explained by the fact that the friction at 
the boundary suppresses the hydrodynamic interaction of the particle with its periodic image. Instead the 
dominant interaction is between the particle and the surface. It is mandatory for a better understanding 
of the system to learn how these finite size effects can be described, quantified, and controlled. 
First we study a system with a constant sphere radius R = 16 and varying system length L. In Fig. [3] 
the drag force F normalized by the Stokes force F$ t = 6n/ivR and the inverse normalized drag force are 
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FIG. 3. Normalized drag force jj- and the inverted normalized drag force versus the normalized distance d/R for different 
system lengths L. The radius of the sphere is fixed at R = 16. The deviation for L = 192 is significant at larger radii, but can be 
neglected for small distances. For L = 512 there is nearly no deviation from the exact solution of Maude \2. It . In addition the 
deviation of the first order approximation 12.2) is below 1% for d > R/2. 



plotted. In the inversed case the deviations for the larger distance d can be seen more clearly. Fig. [3] 
shows that the deviation for the small system L = 192 close to the wall is very small, however the force 
does not converge to the Stokes force F$ t . Here effects similar to the one reported by Hasimoto ( 12.31 > 
appear: for a smaller separation d the force decays with ^ while it approaches a constant value for large 
d. The constant values should be given by the Stokes force F§ t , but can be larger due to the interaction 
with the periodic image. In the l/F plots it can be seen that the deviation is not a constant offset or factor 
but rather starts at a critical value of d/R. From there the force quickly starts to approach a constant 
value. 

In Fig.|4]the relative error E — F ~ F M«'«ie [ s plotted for different system sizes L and a constant radius 
R. The error for the largest system L = 512 in Fig |4]is constantly below 1% for larger distances. At 
distances less than d < R/2 the error rises due to the insufficient resolution of the fluid filled volume 
between the surface of the sphere and the boundary. Another possible effect is the fact that the sphere 
rather jumps over the lattice than performing a continuous movement. Additionally it can be seen that 
for distances less than d = R the error for the different system sizes L collapses. The reason is that for 
smaller distances the lubrication effect which is independent of the system length L dominates the free 
flow and therefore suppresses finite size effects due to the periodic image. The deviation that can be 
seen in the plot at the large d has its origin in the transient. Since the fluid is at rest at the start of the 
simulation and it takes some time to reach a steady state this can only be avoided by longer simulations. 
An interesting fact to point out is that the deviation between the Maude solution and the first order 
approximation is below 1%. 

Since the sphere is discretized on the lattice it is important to understand if this discretization has an 
effect on the drag force. Therefore we perform simulations with a radius of R = 4, 8, 16 at a constant 
ratio R/L = 1 /32 between the radius and the system length. Fig. [5] depicts the normalized drag force 
j- and the inverted normalized drag force ^ versus the normalized distance d/R for different radii R. 
It can be seen that the discretization of the sphere has little influence on the measured force. 

Fig. [6] shows the relative error E for different radii. For all radii the finite size effects due to the 
periodic image are negligible since the ratio between R/L is sufficiently small. The deviation from the 
Maude theory for separations d > R are below 2% for all radii. Therefore one has to concentrate on 
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FIG. 4. Symbols denote the relative error E - 
line shows the first order approximation. As expected the error close to the wall deviates due to the discretization of the small 
distance d. The error becomes larger for large d due to the influence of the periodic image. 




the small distances d where significant deviations appear. In our case this distance is better resolved for 
larger R (note that in the plot the normalized distance is shown). In addition the resolution of the sphere 
is better for larger R. For R = 4 the deviation is more noisy and here the discretization really has an effect 
on the drag force F. However, even a such roughly discretized sphere reproduces the expected result 
surprisingly well. This can probably be explained by averaging out of discretization artefacts due to 
the time and space averaging performed in the simulations. Additionally there should be three or more 
lattice sites between the surface of the sphere and the boundary. If that is not the case the hydrodynamic 
interaction is not resolved sufficiently. If the distance between surface and sphere is smaller than half a 
lattice spacing the two surfaces merge and the method fails. Hence it is advantageous to choose a large 
radius in order to be able to reduce the relative distance to the boundary (in units of the sphere radius) 
or to resolve a possible surface structure. 

For R ^ 8 the deviations have a regular shape and follow the deviation for the first order approxima- 
tion. The trend to follow the first order approximation is stronger for R = 16 but here the noise is reduced 
further and all errors seem to be systematic. Therefore the deviation has to be described as a systematic 
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FIG. 6. Relative error E = F ^Saait ; n % versus the normalized distance d/R for different radii L at a constant ratio R/L = t4t 
(symbols). As expected the error close to the wall deviates, due to the discretization of the small distance d. For R = 4 the 
deviation fluctuates due to the low resolution but for R = 16 it follows the first order approximation. The line corresponds to the 
first order approximation. 



error of the method that has its origin in the "jumping-standing" like movement of the sphere. The first 
order approximation is a quasi static approximation and represents the actual simulated movement more 
correctly than the theory of Maude. 

It should be noted that the finite size effects for R = 16 and R/L — 1/16 due to the interaction 
with the periodic image are much more significant than the discretization effect. By choosing a large 
simulation volume, a radius R > 8 and focusing on the force for separations d < 2R one can reduce 
those effects to a deviation of the measured force from the theoretically predicted value of less than 1 %. 

5. Conclusion 

Lattice Boltzmann simulations of a high-speed drainage of liquid films squeezed between a sphere and 
a smooth no-slip surface have been presented. We have shown that the solution of Maude d2.U can 
be reproduced and demonstrated that at a ratio R/L = 1 /32 finite size effects are below 2% and thus 
can be neglected near the boundary. We have also demonstrated that a sphere radius of R = 8 provides 
a sufficiently well resolved representation of the sphere. Based on this calibration it is possible to 
investigate the influence of different surface properties such as roughness and slip on the drag force on 
AFM based slip measurements. 
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